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Abstract 

Accurate and efficient direct numerical simulation of turbulence in 
the presence of shock waves represents a significant challenge for nu- 
merical methods. The objective of this paper is to evaluate the perfor- 
mance of high order compact and non-compact central spatial differ- 
encing employing total variation diminishing (TVD) shock-capturing 
dissipations as characteristic based filters for two model problems com- 
bining shock wave and shear layer phenomena. A vortex pairing model 
evaluates the ability of the schemes to cope with shear layer instability 
and eddy shock waves, while a shock wave impingement on a spatially- 
evolving mixing layer model studies the accuracy of computation of 
vortices passing through a sequence of shock and expansion waves. 
A drastic increase in accuracy is observed if a suitable artificial com- 
pression formulation is applied to the TVD dissipations. With this 
modification to the filter step the fourth-order non-compact scheme 
shows improved results in comparison to second-order methods, while 
retaining the good shock resolution of the basic TVD scheme. For this 
characteristic based filter approach, however, the benefits of compact 
schemes or schemes with higher than fourth order are not sufficient to 
justify the higher complexity near the boundary and/or the additional 
computational cost. 

‘Queen Mary and Westfield College, London, UK 
fNASA Ames Research Center, Moffett Field, CA, USA 
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1 Introduction 


Time- and space-resolved direct numerical simulation (DNS) of turbulence is 
feasible for many canonical flow problems (Moin & Mahesh, 1998) and with 
increased computer power more applied problems such as separation bubbles 
(Alam & Sandham, 1997) can also be simulated. For many technologically 
important flows compressibility effects are important. Direct simulation is 
already playing an important role in understanding effects such as the re- 
duced growth rate of mixing layers as Mach number is increased (Vreman 
et al., 1996). With high speed compressible flows, however, the potential 
combination of shock waves and turbulent flow represents a significant chal- 
lenge for numerical methods. Although standard total variation diminishing 
(TVD) type of shock-capturing schemes for the Euler equations (Yee, 1989) 
are now routinely used in high speed blast wave simulations with virtually 
non-oscillatory, crisp resolution of discontinuities, for the unaveraged un- 
steady Navier-Stokes equations it was observed (Sandham & Yee 1989) that 
the clipping behavior near extrema of these schemes led to generally poor 
accuracy. 

In response to this difficulty, subsequent large three-dimensional computa- 
tions have either operated at low Reynolds and Mach numbers where the 
shock waves, if present at all, are weak and can be resolved (Luo & Sand- 
ham, 1994) or have used hybrid schemes where shock-capturing schemes are 
turned on only when shock waves are detected (Vreman et al., 1995, Adams 
& Shariff, 1996, Adams, 1997 and Lee et al., 1997). These calculations have 
all used base methods of fourth or sixth-order accuracy, well suited to the 
accurate calculation of turbulence. 

Work has continued on model problems to find improved high order shock- 
capturing schemes. Lumpp (1996a, b) used high order finite volume es- 
sentially non-oscillatory (ENO) schemes for vortex pairing test cases and 
obtained good shock and vortex resolution. However, large grid stencils 
and high computational cost has prohibited application to three-dimensional 
problems. Fu and Ma (1997) have developed a scheme where the group veloc- 
ity near shock waves is fixed by the numerical method such that oscillations 
near shock waves are sucked into the shock wave itself, giving very sharp 
shocks. In another development Yee (1997) has developed a class of compact 
shock-capturing schemes that require smaller grid stencils and operations 
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count than standard schemes. Gustafsson and Olsson (1995) developed sta- 
ble higher-order centered schemes with stable numerical boundary condition 
treatments. For problems containing shocks, Gustafsson and Olsson used 
a scalar shock-capturing filter. Such schemes have advantages over higher- 
order ENO schemes which require very large grid stencils even for modest 
orders of accuracy. (For example, a seven-point grid stencil is required for 
a second-order ENO scheme.) Yee et al. (1998) proposed a combination 
of narrow grid stencil of higher-order classical spatial differencing schemes 
using low order TVD or ENO dissipations as characteristic filters with an 
artificial compression method (ACM) switch. The ACM switch is the same 
as Harten (1978) but applied in a slightly different context. The Yee et al. 
approach is aimed at problems containing vortex convections, shock, shear, 
vortex and turbulence interactions. These characteristic TVD (and ENO) fil- 
ters in conjunction with the artificial compression method can even improve 
fine scale flow structure when applied to existing methods of Yee (1989) and 
Yee (1997). 

A numerical scheme for DNS of shock-turbulence interactions should ideally 
not be significantly more expensive than the classical fourth or sixth-order 
compact or non-compact spatial differencing scheme. It should be possible to 
resolve all scales down to scales of the order of (if necessary) the Kolmogorov 
scales of turbulence accurately and efficiently, while at the same time being 
able to capture steep gradients occurring at smaller scales (e.g. a few mean 
free paths for a strong shock wave). Turbulence mechanisms should not be 
affected by the numerical scheme resulting directly from the governing equa- 
tions. Some early direct simulation codes for incompressible flow were unable 
to sustain turbulence in channel flow due to their properties with respect to 
aliasing errors. Shock-capturing schemes are dissipative and an important 
test of their suitability for turbulence is that they are capable of sustaining 
turbulence. The purpose of the present work is to evaluate a number of can- 
didate schemes discussed in Yee et al. (1998) using model problems. In the 
first of these, a vortex pairing in a time-developing mixing layer, shock waves 
form around the vortices. In the second problem, a shock wave impinging 
on a spatially-evolving mixing layer, the evolving vortices must pass through 
a shock wave, which in turn is deformed by the vortex passage. From these 
two-dimensional tests some conclusions can be drawn regarding numerical 
methods for accurate calculation of vortex motions within an overall shock- 
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capturing framework. 


2 Numerical Techniques 


The numerical methods used for the study are briefly described in this sec- 
tion. The reader is referred to Yee (1989) for a more complete description of 
the basic shock-capturing approaches and Yee et al. (1998) for the complete 
scheme. 

The classical fourth-order Runge-Kutta time discretization and classical fourth 
or sixth-order compact or non-compact spatial differencing as base schemes 
are employed. A second-order TVD dissipation is applied as a filter step at 
the end of the full Runge-Kutta time step in the form of an additional dissi- 
pative numerical flux term 


F j+h ~ 2 Rj+ ^ i+ *' 


( 1 ) 


where R is the right eigenvector matrix of the flux Jacobian from the Euler 
equations and <F is defined by the TVD shock-capturing scheme. The de- 
tailed programming allows the Euler and viscous terms to be computed us- 
ing separate methods. The basic spatial schemes are (i) non-compact central, 
(ii) compact central and (iii) predictor-corrector upwind or upwind biased. 
Non-compact schemes are the standard second, fourth and sixth-order meth- 
ods. Compact schemes are either the standard symmetric fourth-order or 
sixth-order Pade schemes (Ciment & Leventhal, 1975, Hirsh, 1975 and Lele, 
1992). For the purposes of this paper we concentrate on the central schemes. 
Comparable accuracy was obtained with the predictor-corrector upwind or 
upwind biased schemes proposed by Hixon and Turkel (1998). 


2.1 TVD Dissipation as a Filter Step 

The TVD dissipation as a filter step is taken as the diffusive part of an 
upwind explicit shock-capturing scheme of the Harten-Yee type, described in 
Yee (1989). (Other comparable schemes are also applicable.) For a second- 
order upwind TVD scheme the elements of $ i+ i, denoted by are 





( 2 ) 
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The a l .,i are elements of R-}i(Qj+i — Qj), where Q is the vector of conser- 
vative variables, and 

tp(z) = ^/(6 + z 2 ), (3) 

and , 

, _ (4l 

7 "* 2 (<■$♦»>’+« 

To avoid an additional logical statement in the actual coding, e is added to 
the 7* . The a 1 , i are characteristic speeds of the Euler flux evaluated using 

Roe’s average. In all of the computations, we take e = 10 -7 . We use the flux 
limiter given by 



a LiI(«i+i) 2 + e l + + e l 


i+5 


(«' + i) 2 + K_i) +2e 


( 5 ) 


Elements of R j+ i are computed using Roe’s averaging procedure. In all of 
the computations, the value of S was taken to be 1/16 to satisfy an entropy 
condition. The resolution of the fine scale flow structure showed minor sen- 
sitivity to the value of this constant. See Yee et al. (1998) for a discussion. 

2.2 The ACM Switch 

An artificial compression modification can be made to the TVD dissipation, 
hereafter, referred to as “the ACM/TVD filter”. The term (/>j+ 1/2 mu hi~ 
plied by a switch designed by Harten (1978) 


icmax(0*, 0 l j+1 ), (6) 


where 




l<X- 


■3+ 1/2 


- a' 


I-1/2! 


. (7) 

K+1/2I + K-1/2I + € 

This serves to limit the action of the TVD dissipation to the immediate vicin- 


ity (within one grid cell) of the discontinuity, as detected by steep gradients 
in the characteristic variables. The constant e = 10 -7 is again used to avoid 
any potential singularity without an extra logical statement in the actual 
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Method 

Order 

(Euler) 

Order 

(viscous) 

Shock- 

capturing 

Artificial 

compression 

Compact 

CEN44 

4 

4 

No 

No 

No 

TVD22 

2 

2 

Yes 

No 

No 

TVD44 

4 

4 

Yes 

No 

No 

TVD66 

6 

6 

Yes 

No 

No 

ACM22 

2 

2 

Yes 

Yes 

No 

ACM44 

4 

4 

Yes 

Yes 

No 

ACM66 

6 

6 

Yes 

Yes 

No 

ACM44C 

4 

4 

Yes 

Yes 

Yes 

ACM66C 

6 

6 

Yes 

Yes 

Yes 


Table 1: Notation for numerical methods. Order of accuracy refers to the 
formal order of the base scheme. 


coding. The optimum value of k is problem-dependent and will be discussed 
further in a later section. See Yee et al. (1998) for additional discussion. 

The various combinations of schemes that will be used in this paper are 
shown on Table 1. The notation “TVD” with the various orders attached at 
the end means the second-order TVD scheme dissipation (without the ACM 
switch) is used as the filter with the indicated order of the base scheme for the 
convection and viscous terms. For simplicity of discussion, unless otherwise 
indicated, the term TVD or ACM scheme means the selected base schemes 
indicated in Table 1 using the TVD or ACM/TVD filter. 


3 Test Cases and Reference Solutions 

The performance of the candidate schemes is based on a Cartesian grid with 
an analytical stretching in the y-direction only. For the present study no 
attempt is made to use a more sophisticated adaptive grid. 

3.1 Vortex Pairing in a Time-Developing Mixing Layer 

The first test case is of vortex growth and pairing in a temporal mixing layer 
at a convective Mach number equal to 0.8. At this Mach number there are 
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shock waves (shocklets) that form around the vortices and the problem is to 
compute accurately the vortex evolution while avoiding oscillations around 
the shocks. Previous calculations of the problem can be found in Sandham 
and Reynolds (1989), Lumpp(1996a,b) and Fu and Ma (1997). Here we set 
up a base flow as in Sandham and Yee (1989) 

u = 0.5 tanh(2y), (8) 


with velocities normalized by the velocity jump Ui — U 2 across the shear layer 
and distances normalized by vorticity thickness, 


. U 1 -U 2 

( du / dy ) max 


( 9 ) 


Subscripts 1 and 2 refer to the upper ( y > 0) and lower (y < 0) streams of 
fluid respectively. The normalized temperature and hence local sound speed 
squared is determined from an assumption of constant stagnation enthalpy 


c 2 = cl + 2— Vf - u 2 ). 


( 10 ) 


Equal pressure through the mixing layer is assumed. Therefore, for this 
configuration of U 2 — — U \ both fluid streams have the same density and 
temperature for y — > ioo. The Reynolds number defined by the velocity 
jump, vorticity thickness and kinematic viscosity at the free-stream temper- 
ature is set here to be 1000. The Prandtl number is set to 0.72, the ratio of 
specific heats is taken as 7 = 1.4 and Sutherland’s law 

_H_ (c 2 /4)i- 5 (l + 110.3 /Tr) 

Hr c 2 /c 2 r + 110.3/T* 

with reference temperature T R = 300 K is used for the viscosity variation 
with temperature. The reference sound speed squared c 2 R is taken as the 
average of c 2 over the two free streams with hr the viscosity corresponding 
to c 2 r . 

Disturbances are added to the velocity components in the form of simple 
waves. For the normal component of velocity we have the perturbation 


2 

v' = a k cos(2n kx/L x + <f> k ) exp (~y 2 /b), 

k = 1 


( 12 ) 
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where L x = 30 is the box length in the a;— direction and b = 10 is the 
y— modulation. In our test case we simulate pairing in the center of the 
computational box, by choosing the initially most unstable wave k = 2 to 
have amplitude a 2 = 0.05 and phase fa = -vr/2, and the subharmonic wave 
k = 1 with oi = 0.01 and fa = -tt/ 2. The u-velocity perturbations are 
found by assuming that the total perturbation is divergence free. These 
fluctuations correspond only approximately to eigenfunctions of the linear 
stability problem for a compressible mixing layer, but they serve the purpose 
of initiating the instability of the mixing layer and have the advantage as a 
test case in that they can be easily coded. 

Numerically the grid is equally spaced and periodic in the x direction and 
stretched in the y-direction, using the mapping 

_ Ly smh(b y ri) , 13 \ 

^ 2 sinh(6 y ) 

where we take the box size in the y — direction Ly — 100, and the stretching 
factor by = 3.4. The mapped coordinate rj is equally spaced and runs from 
-1 to +1. The boundaries at ±L y / 2 are taken to be slip walls. For example, 
at the lower boundary 


Pi 

= P2, 

(14) 

(pu) 1 

- (pu) 2 , 

(15) 

(pv) 1 

.= 0 , 

(16) 

(Et) i 

= [4(£t) 2 — (£t) 3 ]/3, 

(17) 


where subscripts here refer to the grid point and E T is the total energy. 

Sample results for this configuration were obtained from computations on 
a 201 x 201 grid, using the fourth-order non-compact central spatial differ- 
encing, and the ACM/TVD filter (ACM44). The constant /c=0.7 for the 
u ± c and v ± c nonlinear characteristic fields and /c=0.35 for the u and v 
linear characteristic fields were used for the computations. Figure la shows 
snapshots of the pressure field at times t = 40, 80, 120, 160, illustrating the 
roll-up of the primary vortices, followed by vortex merging. Figure lb - Id 
show the corresponding density, vorticity and temperature contours. Shock 
waves form around the vortices, with a peak Mach number ahead of the vor- 
tex of approximately 1.55 at t = 120. Figure 2 shows the growth of two 


8 



measures of shear layer thickness, the vorticity thickness defined by equation 
(11) and a momentum thickness, defined here by 

B = [ Ly/2 p(Ui -u)(u- U 2 )dy, (18) 

J ~I J y /2 

with averages being taken over the periodic direction x. Both thickness mea- 
surements show rapid growth up to a peak at around t = 120 followed by 
decay. The vorticity thickness is generally more sensitive to details of the flow 
evolution. Because of the finite box in x the flow would eventually relami- 
narize. This case is used as a reference high resolution case for comparison of 
schemes in Section 4. It was checked to be adequately converged by running 
another simulation with 401 x 401 grid points. 

3.2 Shock Wave Impingement on a Spatially-Evolving 
Mixing Layer 

The second test case has been developed to test the behavior of the schemes 
for shock waves interacting with shear layers where the vortices arising from 
shear layer instability are forced to pass through a shock wave. An oblique 
shock is made to impact on a spatially-developing mixing layer at an initial 
convective Mach number of 0.6. The shear layer vortices pass through the 
shock system and later through another shock, imposed by reflection from a 
(slip) wall at the lower boundary. The problem has been arranged so that 
the Mach number at the outflow boundary is everywhere supersonic, so no 
explicit outflow boundary conditions are required. This allows us to focus on 
properties of the numerical schemes rather than on the boundary treatment. 

Figure 3 illustrates the nature of the flow on a 321 x 81 grid. This result was 
taken from a laminar flow simulation with no incoming perturbations. The 
shear layer originates at x — 0, y = 0 at the center of the left hand boundary. 
An oblique shock originates from the top left hand corner and this impacts 
on the shear layer at around x = 90. The shear layer is deflected by the 
interaction. Afterwards we have a shock wave below the shear layer and an 
expansion fan above it. The shock wave reflects from the lower solid wall and 
passes back through the shear layer. The lower wall uses a slip condition so 
no viscous boundary layer forms and we focus on the shock- wave interaction 
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Property 

(1) 

(2) 

( 3 ) 

( 4 ) 

( 5 ) 

u — velocity 

3.0000 

2.0000 

2.9709 

2.9792 

1.9001 

velocity 

0.0000 

0.0000 

-0.1367 

-0.1996 

-0.1273 

9 (degrees) 

0.0000 

0.0000 

2.6343 

3.8330 

3.8330 

density p 

1.6374 

0.3626 

2.1101 

1.8823 

0.4173 

pressure p 

0.3327 

0.3327 

0.4754 

0.4051 

0.4051 

sound speed c 

0.5333 

1.1333 

0.5616 

0.5489 

1.1658 

Mach number \M\ 

5.6250 

1.7647 

5.2956 

5.4396 

1.6335 


Table 2: Flow properties for the shock- wave/shear-layer test case in various 
regions of the flow: (1) upper stream inflow, (2) lower stream inflow, (3) 
upper stream after oblique shock, (4) upper stream after expansion fan, (5) 
lower stream after shock wave. 


with the unstable shear layer. The full no-slip problem would, however, make 
a challenging test case for the future. 

The inflow is specified again with a hyperbolic tangent profile, this time as 

u = 2.5 -I- 0.5 tanh(2y), (19) 

giving a mixing layer with upper velocity U\ = 3, lower velocity C /2 = 2, and 
hence a velocity ratio of 1.5. Equal pressures and stagnation enthalpies are 
assumed for the two streams, with convective Mach number, defined by 

M c - U '~ U \ (20) 

Ci 4- C 2 

where ci and C 2 are the free stream sound speeds, equal to 0.6. The reference 
density is taken as the average of the two free streams and a reference pressure 
as {pi + p 2 )(C/i — C/2) 2 /2. This allows one to compute the inflow parameters 
as given in the first two columns of Table 2. Inflow sound speed squared is 
found from the relation for constant stagnation enthalpy (12). 

The upper boundary condition is taken from the flow properties behind an 
oblique shock with angle /? = 12°, given in column 3 of Table 2. The table 
also gives the properties behind the expansion fan (column 4) and after the 
oblique shock on the lower stream of fluid (column 5) computed by standard 
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gasdynamics methods with (3 = 38.118°. The conditions in regions 4 and 
5 do not correspond exactly to the simulations due to the finite thickness 
of the shear layer in practice. The Mach number of the lower stream after 
this shock is approximately M 5 = 1.6335 and remains supersonic through all 
the successive shocks and expansion fans up to the outflow boundary. The 
resulting shock waves are not strong, but tests showed that they could not 
be computed without using shock-capturing techniques. The lower boundary 
was specified with the same slip condition used for the pairing case (Equations 
(16-19)). 

The Prandtl number and ratio of specific heats were taken to be the same as 
for the vortex pairing test case. The Reynolds number was chosen to be 500. 

Fluctuations are added to the inflow as 

v' = 53 a,k cos(2n kt/T + (f>k) exp (—y 2 /b), (21) 

k= 1 

with period T = \/U c , wavelength A = 30, convective velocity U c = 2.68 
(defined by U c = ( U\C 2 + U 2 C\)/{c\ + c 2 )) an d b = 10. For k = 1 we take 
ai = 0.05 and 0 = 0, and for k = 2 we take a 2 = 0.05 and 0 = 7r/2. No 
perturbations are added to the u — component of velocity. 

The grid is taken to be uniform in x and stretched in y according to equation 
(15) with by = 1. This stretching is much milder than for the pairing problem, 
as we have to resolve the shear layer even when it deflects away from y — 0. 
The box lengths were taken to be L x = 200 and L y = 40. 

A reference solution for the pressure, density and temperature fields is shown 
on figure 4 for a computation on a grid of 641 x 161. The computation employs 
the fourth-order central differencing as the base scheme in conjunction with 
the ACM/TVD filter (ACM44) with k = 0.35 for nonlinear characteristic 
fields and k = 0.175 for the linear characteristic fields. The vortex cores are 
located by low pressure regions and the stagnation zones between vortices 
by high pressure regions. The shock waves are seen to be deformed by the 
passage of the vortices. Another interesting observation is the way the core 
of the vortex at x = 148 has been split into two by its passage through the 
reflected shock wave. In spite of the relatively high amplitude chosen for the 
subharmonic inflow perturbation there is no pairing of vortices within the 
computational box. We do, however, begin to see eddy shock waves around 


11 



the vortices near the end of the computational box where the local convective 
Mach number has increased to around 0.66. The oscillations seen near the 
upper boundary for x > 120 occur where the small Mach waves from the 
initial perturbations arrive at the upper boundary. The use of characteristic 
boundary conditions should remove this problem. Practically, the amplitude 
of oscillations is not sufficient to cause numerical instability or affect the 
remainder of the flow. 


4 Computational Results 

The notation that will be used for discussing the results for different numer- 
ical schemes is shown on Table 1. To examine the resolution of the proposed 
schemes where shock waves are absent, the computation is compared with 
the CEN44 (the classical spatially fourth-order central differencing for the 
convection and diffusion terms) before shock waves were developed for the 
vortex pairing case. Good agreements were obtained. 

The performance of these schemes with the presence of shock waves and 
turbulence is evaluated based on the following factors: 

(a) Effect of the ACM term 

(b) Effect of the order of the base scheme 

(c) Effect of the grid size (grid refinement study) 

(d) Effect of employing a compact or non-compact base scheme 

(e) Effect of the adjustable constant k for the particular physics 

(f) Shear and fine flow structure capturing capability 

4.1 Vortex Pairing 

For an initial comparison of the schemes we compute the test case of Section 

3.1 on a grid of 101 x 101 with a time step A t — 0.1, running up to t — 160. 

Effect of the A CM Term: 

Figure 5 shows the effect of the artificial compression method for the vortic- 
ity and momentum thickness variation with time for the TVD44 and ACM44 
schemes compared with the reference results from a fine grid (Section 3.1). 
Results from the ACM method are far superior to those from the standard 
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TVD formulation. Note that there is no improvement in the shock resolu- 
tion since the ACM term limited the amount of dissipation away from high 
gradient areas whereas the shock resolution is dictated by the flux limiter. 

Effect of the Order of the Base Scheme : 

Figure 6 shows the effect of increasing accuracy from second to fourth and 
sixth order using the TVD filter (TVD22, TVD44 and TVD66). As can be 
seen there is almost no improvement as the order of accuracy is raised. Fig- 
ure 7 shows the same plot for the ACM/TVD filter (ACM22, ACM44 and 
ACM66). Here there is an improvement, although the results even for the 
lowest order are quite good. There is little to choose in the shock resolution 
properties with the variation in order of accuracy. We choose to compare 
temperature contours, which are most sensitive to oscillations (Lumpp, pri- 
vate communication). Figure 8 shows a comparison of the ACM schemes of 
various order with the reference solution. All the schemes here capture the 
shock waves with minimal oscillations. The temperature contours for the 
TVD filter of the various order using a 101 x 101 grid, although not shown, 
are not even nearly as accurate as the ACM44 using a 41 x 41 grid. See the 
last plot of Figure 15 for a comparison. It can be seen that there is a sig- 
nificant advantage in moving from second to fourth order, but less is gained 
in moving from fourth to sixth order using TVD or ACM/TVD as a filter. 
This is in contrast to an isentropic vortex convection test case where there 
are definite benefits of moving from fourth order to sixth order (Yee et al., 
1998). It appears that the effect of order of accuracy are more pronounced 
for long time integrations of pure convection. 

Effect of the Grid Size ( Grid Refinement Study): 

To investigate the effect of order of the accuracy in more detail we consider 
simulations on a very coarse grid of 41 x 41 points. Such a case corresponds 
in practice to simulation of scales of turbulence arising from shear layers only 
two or three computational cells across. Figure 9 shows results for the ACM 
schemes. To ensure that the fine scale flow structure is fully resolved by the 
reference grid 201 x 201, the same simulation was done on a 401 x 401 grid 
(figures not shown). 
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Effect of Compact or Non- Compact Base Scheme : 

For wave propagation and computational problems the performance of fourth 
and sixth-order compact schemes, although more CPU intensive, appears to 
be superior to their non-compact cousin. For problems with shock waves the 
benefit of compact over non-compact schemes is less known due to the filter 
step. Figure 10 shows results for the fourth and sixth-order compact schemes, 
which are similar to results from the sixth-order non-compact scheme. Again 
there is little improvement compared with the fourth-order non-compact 
scheme. A conclusion is that the use of the ACM/TVD in the filter step 
is essential to get the benefits of moving from second to fourth order, but 
even with this method there is little benefit in moving to even higher-order 
schemes. 

Effect of the Adjustable Constant k : 

The ACM switch of Harten has been demonstrated to give good shock resolu- 
tion and to be essential if the benefits of higher-order discretization schemes 
are to be realized. There is, however, an adjustable parameter k in the for- 
mulation, and results are sensitive to the precise choice of its value. Figure 
11 illustrates the effect on the result for the pairing test case of reducing the 
parameter from 0.7 to 0.35. The vorticity and momentum thickness devel- 
opment is improved due to the reduction in numerical dissipation. From the 
temperature contours on Figure 12 it can be seen that this has been achieved 
at the cost of formation of small oscillations around the shock wave. For the 
present problem one would be ready to pay this price to get the more accurate 
vortex evolution. However, in general it is not known how such numerically- 
induced oscillations interact with small scales of turbulence. For the current 
method the correct procedure for a simulation of shock-turbulence interac- 
tion would be to find the smallest value of k to resolve the shock waves 
satisfactorily and then increase the grid resolution until the turbulence is ad- 
equately resolved. There are perhaps other formulations of the ACM switch 
parameter k based on the flow physics that can perform the adjustment to 
higher values automatically when stronger shock waves are present. This is 
a subject of future research. 

Shear and Fine Flow Structure Capturing Capability : 

Shock-capturing schemes are designed to accurately capture shock waves, but 
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with a less accurate capturing capability for contact discontinuities. In fact, 
the mixing layer seen at a large scale is a contact discontinuity. If one uses 
enough grid points to resolve the region of high shear in conjunction with 
physical viscosity, it might not need to be ‘captured’. Contact discontinuities 
relate to the characteristic velocities u and v. As an experiment these linear 
characteristics were computed with no numerical dissipation. Filters were 
only applied to the nonlinear characteristics fields u± c and v ± c using the 
ACM44 method. Interestingly, the computation was no less stable than that 
the full TVD or ACM schemes (applying k to all of the characteristic fields). 
Figure 13 shows the vorticity thickness evolution and Figure 14 the resulting 
temperature field. It can be seen that good results were obtained, although 
there is a trace of oscillation near the shock wave. This could be remedied 
by increasing k slightly. See Yee et al. (1998) for an illustration. However, 
the flow features of the shear and fine flow structure are accurately captured 
with similar resolution as the 201 x 201 grid with ACM44 applied to all of 
the characteristic fields (Figure 15a). 

To balance the shear and shock-capturing, one alternative is to switch to a 
more compressive limiter (see Yee (1989)) for the linear characteristics fields. 
Another alternative is to reduce the value of k for the u, v linear fields. The 
comparison of using different values of k for the linear and nonlinear fields is 
also shown in Figures 13 and 14. 

Figure 15 shows the comparison of the different k values for the three grids. 
Figure 16 shows the comparison of the five classical flux limiters (see Equa- 
tions (4.34c)-(4.34g) of Yee (1989)) using ACM44. 

4.2 Shock-Shear Layer Interaction 

The test case defined in Section 3.2 was run on a grid of 321 x 81 with 
At = 0.12 up to t = 120 with k = 0.35. Figure 17 shows the density field 
for the TVD44 and ACM44 schemes compared with the reference solution 
(641 x 161 grid). Ag ain it can be seen that the ACM modification is essential 
for obtaining good vortex evolution (additionally better shock resolution is 
obtained). For a more quantitative comparison Figures 18 and 19 show slices 
through the density field between x = 50 and x = 150 at y = 0 (the centerline 
of the computational domain). From Figure 18 it is apparent that all the 
standard TVD schemes, of whatever order, miss the correct vortex formation. 
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Method 

Method 

Pairing 
(200 steps) 

Shock-wave 
(100 steps) 

CEN44 

19.5 

22.1 

TVD22 

18.8 

22.8 

TVD44 

23.4 

26.9 

ACM22 

19.6 

23.8 

ACM44 

24.2 

27.9 


Table 3: CPU times for various numerical options, measured in Cray C90 
seconds. For a key to the various methods see Section 4.1. 


From Figure 19 there are some visible benefits in moving from second to 
higher-order differencing, both in the amplitude of the fluctuations and in 
the correct convective velocity of the vortices. Figures 20 and 21 compare 
the density contour using limiter 3 and limiter 5. These figures also compare 
the use of different k values for the linear and nonlinear characteristic fields. 

4.3 Computational Costs 

The research code that we built consists of many schemes and options and 
thus is far from being an optimum code. Based on this fact, the CPU times 
for the Cray C90 are given on Table 3 for the two model problems. Times 
are given in seconds for 200 steps of the vortex pairing test case on a grid 
of 101 x 101 and for 100 steps of the shock-wave impingement problem on 
a grid of 321 x 81. 

The non-compact base schemes with the ACM/TVD filter are only around 
25% more expensive than the same base schemes without ACM/TVD filter. 
This has been achieved by only requiring one application of the ACM/TVD 
filter per full time step for the convection terms. 

For the Cray C90 it was found that the compact schemes were significantly 
more expensive, but this result is distorted for the present code by incomplete 
vectorization. An extra cost of around 2/3 is expected from considerations 
of operations count. 
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5 Conclusions 


Two model problems were considered to test the behavior of shock-capturing 
schemes in flows where the evolution of vortices from shear layer instabilities 
plays a crucial role. Standard total variation diminishing schemes do not 
perform well for these problems. It was found that an artificial compression 
modification gave much improved results and with this approach there are 
definite benefits in moving from second to fourth-order differencing schemes. 
The benefits were less pronounced in moving to sixth order. The fourth order 
ACM44 schemes show highly accurate shock-wave and shear-layer capability. 
The non-compact fourth-order base scheme requires a five point stencil, the 
same as the second-order TVD scheme. There is potentially an additional 
problem near boundaries in that reduced order and/or upwind schemes must 
be used. Stable boundary conditions such as those proposed by Gustafs- 
son and Olsson (1995) must be applied. If an adaptive grid is used, higher 
accuracy can be obtained. 
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Figure 7. The effect of order of accuracy on the ACM methods, ACM22 
ACM44 ( — ) and ACM66 (chain dotted), compared with the reference 
solution (solid line): (a) vorticity thickness 5 U , (b) momentum thickness 9 
for a 101 x 101 grid. 
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